options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -4384.8 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24      -36.8162   0.000474379    0.00124184      0.9102      0.9102       43    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -36.82
##    b0         6.43
##    b1         1.75
##    s          3.82
##    m0[1]     28.17
##    m0[2]     35.44
##    m0[3]     10.55
##    m0[4]     11.11
##    m0[5]     14.44
##    m0[6]     20.73
##    m0[7]     43.22
##    m0[8]      4.41
##    m0[9]     35.08
##    m0[10]    31.91
##    m0[11]    37.97
##    m0[12]    31.78
##    m0[13]    29.86
##    m0[14]    39.48
##    m0[15]    30.89
##    m0[16]    45.01
##    m0[17]    16.60
##    m0[18]     9.70
##    m0[19]     4.38
##    m0[20]    25.93
##    y0[1]     23.02
##    y0[2]     31.09
##    y0[3]     14.82
##    y0[4]      9.54
##    y0[5]     15.89
##    y0[6]     23.38
##    y0[7]     41.94
##    y0[8]      5.91
##    y0[9]     33.80
##    y0[10]    23.17
##    y0[11]    42.86
##    y0[12]    36.30
##    y0[13]    32.92
##    y0[14]    35.46
##    y0[15]    36.00
##    y0[16]    43.62
##    y0[17]    23.45
##    y0[18]    12.53
##    y0[19]     7.67
##    y0[20]    22.98
##    m1[1]      4.38
##    m1[2]      8.90
##    m1[3]     13.41
##    m1[4]     17.92
##    m1[5]     22.44
##    m1[6]     26.95
##    m1[7]     31.47
##    m1[8]     35.98
##    m1[9]     40.49
##    m1[10]    45.01
##    y1[1]      0.92
##    y1[2]     12.45
##    y1[3]      9.97
##    y1[4]     19.15
##    y1[5]     25.82
##    y1[6]     20.50
##    y1[7]     32.22
##    y1[8]     29.03
##    y1[9]     43.35
##    y1[10]    41.12
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -36.99 -36.71 1.21 0.99 -39.39 -35.65 1.00      667      588
##    b0       6.40   6.44 1.66 1.61   3.69   8.95 1.00      328      345
##    b1       1.75   1.75 0.13 0.13   1.54   1.96 1.00      449      573
##    s        4.31   4.24 0.76 0.75   3.26   5.59 1.01     1081     1294
##    m0[1]   28.14  28.16 1.00 0.95  26.47  29.79 1.00     1544     1293
##    m0[2]   35.40  35.42 1.25 1.21  33.33  37.41 1.00     2487     1410
##    m0[3]   10.52  10.53 1.42 1.38   8.22  12.70 1.00      333      374
##    m0[4]   11.08  11.10 1.39 1.36   8.84  13.24 1.00      335      391
##    m0[5]   14.40  14.45 1.22 1.20  12.43  16.33 1.00      354      433
##    m0[6]   20.70  20.71 1.00 0.99  19.07  22.35 1.00      494      578
##    m0[7]   43.18  43.21 1.68 1.65  40.40  45.83 1.00     1376     1453
##    m0[8]    4.38   4.43 1.78 1.77   1.46   7.13 1.00      328      348
##    m0[9]   35.04  35.07 1.24 1.19  32.99  37.03 1.00     2524     1410
##    m0[10]  31.87  31.89 1.10 1.08  30.04  33.68 1.00     2510     1526
##    m0[11]  37.94  37.95 1.38 1.35  35.62  40.07 1.00     2055     1506
##    m0[12]  31.75  31.76 1.10 1.07  29.93  33.55 1.00     2495     1526
##    m0[13]  29.82  29.84 1.04 0.99  28.14  31.52 1.00     2069     1354
##    m0[14]  39.44  39.46 1.46 1.43  36.99  41.72 1.00     1812     1369
##    m0[15]  30.86  30.86 1.07 1.02  29.10  32.61 1.00     2346     1453
##    m0[16]  44.97  45.01 1.80 1.75  42.01  47.79 1.00     1240     1342
##    m0[17]  16.57  16.62 1.13 1.09  14.75  18.38 1.00      379      486
##    m0[18]   9.67   9.69 1.46 1.42   7.28  11.93 1.00      331      364
##    m0[19]   4.35   4.40 1.78 1.77   1.43   7.10 1.00      328      348
##    m0[20]  25.89  25.90 0.96 0.93  24.31  27.51 1.00     1017     1183
##    y0[1]   28.34  28.41 4.52 4.33  21.00  35.70 1.00     1819     1919
##    y0[2]   35.33  35.46 4.70 4.55  27.18  43.14 1.00     2070     1765
##    y0[3]   10.57  10.67 4.59 4.44   2.89  18.02 1.00     1506     1999
##    y0[4]   10.95  10.95 4.56 4.25   3.35  18.43 1.00     1485     1643
##    y0[5]   14.32  14.38 4.49 4.37   6.66  21.40 1.00     1445     1513
##    y0[6]   20.78  20.83 4.38 4.26  13.42  27.86 1.00     2036     1972
##    y0[7]   43.14  42.98 4.73 4.66  35.51  50.86 1.00     1825     1553
##    y0[8]    4.49   4.51 4.83 4.71  -3.46  12.40 1.00     1247     1542
##    y0[9]   34.89  34.81 4.49 4.33  27.56  42.38 1.00     2041     1804
##    y0[10]  31.75  31.78 4.58 4.20  24.20  39.25 1.00     1966     2011
##    y0[11]  37.97  37.98 4.51 4.33  30.69  45.34 1.00     2046     1856
##    y0[12]  31.74  31.71 4.44 4.25  24.39  39.16 1.00     1857     1914
##    y0[13]  29.79  29.82 4.45 4.32  22.35  37.12 1.00     2063     2056
##    y0[14]  39.47  39.33 4.54 4.53  32.11  47.21 1.00     2055     1970
##    y0[15]  30.78  30.81 4.48 4.38  23.36  38.00 1.00     2038     1819
##    y0[16]  44.89  44.99 4.84 4.70  36.83  52.78 1.00     1866     1932
##    y0[17]  16.55  16.49 4.55 4.34   8.99  23.82 1.00     1640     1918
##    y0[18]   9.87   9.92 4.57 4.45   2.24  17.18 1.00     1341     1442
##    y0[19]   4.45   4.46 4.67 4.63  -3.13  12.08 1.00     1133     1427
##    y0[20]  25.98  25.95 4.55 4.35  18.50  33.57 1.00     2042     1927
##    m1[1]    4.35   4.40 1.78 1.77   1.43   7.10 1.00      328      348
##    m1[2]    8.87   8.90 1.51 1.47   6.40  11.21 1.00      330      369
##    m1[3]   13.38  13.41 1.27 1.23  11.32  15.37 1.01      346      411
##    m1[4]   17.89  17.92 1.08 1.06  16.15  19.62 1.00      404      537
##    m1[5]   22.41  22.42 0.97 0.95  20.83  24.01 1.00      595      748
##    m1[6]   26.92  26.93 0.98 0.94  25.31  28.55 1.00     1230     1196
##    m1[7]   31.43  31.44 1.09 1.05  29.64  33.23 1.00     2447     1525
##    m1[8]   35.94  35.97 1.28 1.24  33.82  37.98 1.00     2535     1393
##    m1[9]   40.46  40.48 1.52 1.49  37.92  42.84 1.00     1672     1289
##    m1[10]  44.97  45.01 1.80 1.75  42.01  47.79 1.00     1240     1342
##    y1[1]    4.48   4.47 4.62 4.42  -2.85  12.06 1.00     1057     1380
##    y1[2]    8.91   8.95 4.57 4.31   1.38  16.46 1.00     1584     1955
##    y1[3]   13.43  13.34 4.59 4.37   6.12  20.94 1.00     1703     1876
##    y1[4]   17.78  17.61 4.59 4.34  10.36  25.01 1.00     1565     1899
##    y1[5]   22.57  22.61 4.54 4.33  15.06  29.91 1.00     1664     1855
##    y1[6]   26.88  27.01 4.55 4.19  19.33  34.40 1.00     2017     2001
##    y1[7]   31.36  31.41 4.55 4.27  24.03  38.85 1.00     2017     2083
##    y1[8]   35.81  35.81 4.53 4.55  28.39  43.14 1.00     1943     1678
##    y1[9]   40.44  40.38 4.66 4.52  32.78  48.02 1.00     1814     1917
##    y1[10]  44.91  44.80 4.68 4.63  37.44  52.38 1.00     2113     1974

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -41.21 -43.66 12.40 10.50 -57.54 -16.72 1.11       32       61
##    b0       6.51   6.48  1.81  1.70   3.58   9.54 1.01     1065     1471
##    b1       1.74   1.75  0.14  0.13   1.51   1.97 1.01      903     1072
##    s        3.37   3.46  1.30  1.28   1.17   5.39 1.03      173      239
##    sx       1.39   1.29  0.91  1.06   0.20   2.86 1.09       36       75
##    x0[1]   13.98  13.49  1.76  1.67  11.90  17.09 1.03      179     1011
##    x0[2]   16.54  16.55  1.06  0.75  14.77  18.28 1.02     2702     1516
##    x0[3]    1.56   1.87  1.29  0.99  -0.88   3.19 1.01      458     1100
##    x0[4]    3.81   3.46  1.46  1.35   1.98   6.46 1.02      174     1440
##    x0[5]    4.52   4.54  1.06  0.76   2.77   6.24 1.04     2648     1775
##    x0[6]    7.25   7.56  1.29  1.13   4.87   8.89 1.02      250     1150
##    x0[7]   21.31  21.16  1.18  0.83  19.64  23.35 1.03     1465      948
##    x0[8]   -0.79  -0.93  1.16  0.86  -2.61   1.21 1.02      769     1578
##    x0[9]   16.59  16.49  1.07  0.76  14.97  18.35 1.03     2167     1335
##    x0[10]  13.70  13.95  1.26  1.14  11.41  15.35 1.02      364     1903
##    x0[11]  18.80  18.51  1.31  1.04  17.17  21.13 1.01      466      973
##    x0[12]  14.72  14.63  1.06  0.79  13.05  16.57 1.03     1784     1287
##    x0[13]  13.70  13.60  1.07  0.83  12.05  15.52 1.02     1501     1435
##    x0[14]  19.62  19.33  1.31  0.98  17.89  21.93 1.02      486     1111
##    x0[15]  12.70  13.09  1.52  1.45   9.89  14.55 1.03      139     1158
##    x0[16]  20.96  21.32  1.44  1.29  18.38  22.74 1.02      185     1288
##    x0[17]   4.90   5.22  1.33  1.12   2.46   6.59 1.02      358     1201
##    x0[18]   0.90   1.27  1.38  1.12  -1.70   2.55 1.02      365      983
##    x0[19]  -0.76  -0.92  1.14  0.89  -2.48   1.23 1.03      764     1527
##    x0[20]  12.18  11.84  1.41  1.25  10.41  14.65 1.02      276     1413
##    m0[1]   30.80  30.28  2.99  3.23  26.67  35.70 1.02      208     1579
##    m0[2]   35.28  35.29  1.91  1.68  31.99  38.32 1.01     3888     2506
##    m0[3]    9.28   9.38  2.27  2.32   5.62  12.90 1.02      343     1969
##    m0[4]   13.16  12.89  2.64  2.82   9.22  17.58 1.02      235     1622
##    m0[5]   14.40  14.40  1.91  1.67  11.30  17.51 1.01     3094     2649
##    m0[6]   19.16  19.40  2.28  2.39  15.29  22.57 1.02      170     1706
##    m0[7]   43.54  43.64  2.07  1.91  39.99  46.88 1.01     1784     2614
##    m0[8]    5.19   5.23  2.20  2.08   1.57   8.73 1.00     1141     2147
##    m0[9]   35.35  35.37  1.89  1.68  32.25  38.40 1.01     3345     2282
##    m0[10]  30.34  30.51  2.27  2.41  26.55  33.70 1.02      276     2416
##    m0[11]  39.19  39.10  2.26  2.31  35.67  42.86 1.01      512     2526
##    m0[12]  32.12  32.13  1.88  1.60  28.97  35.16 1.01     2432     2278
##    m0[13]  30.34  30.34  1.86  1.72  27.29  33.36 1.01     1756     2645
##    m0[14]  40.60  40.58  2.22  2.26  36.96  44.15 1.01      518     2563
##    m0[15]  28.61  29.11  2.71  2.93  24.01  32.39 1.03      105     1631
##    m0[16]  42.95  43.09  2.66  2.99  38.65  47.03 1.02      182     2215
##    m0[17]  15.08  15.27  2.34  2.53  11.23  18.63 1.02      310     2050
##    m0[18]   8.13   8.26  2.39  2.56   4.27  12.00 1.02      338     1604
##    m0[19]   5.24   5.26  2.18  2.07   1.65   8.65 1.01      865     1767
##    m0[20]  27.69  27.36  2.43  2.55  24.11  31.77 1.02      303     2359
##    y0[1]   30.75  31.15  4.70  4.76  22.72  37.48 1.01      425     2952
##    y0[2]   35.33  35.30  4.11  3.69  28.53  42.22 1.00     4227     3618
##    y0[3]    9.32   8.91  4.36  3.96   2.69  16.86 1.00     1659     2918
##    y0[4]   13.10  13.50  4.48  4.34   5.27  19.71 1.01      693     2926
##    y0[5]   14.41  14.34  4.16  3.53   7.71  21.36 1.00     3622     2945
##    y0[6]   19.14  18.78  4.20  3.91  12.87  26.36 1.01      959     3038
##    y0[7]   43.55  43.64  4.09  3.48  36.85  50.14 1.00     3382     2902
##    y0[8]    5.14   5.36  4.20  3.75  -2.16  11.63 1.00     2523     2828
##    y0[9]   35.29  35.38  4.03  3.53  28.42  41.86 1.00     3708     3299
##    y0[10]  30.23  29.88  4.28  3.90  23.91  37.67 1.01     1501     3525
##    y0[11]  39.18  39.49  4.18  3.76  31.81  45.52 1.00     1763     2431
##    y0[12]  32.18  32.36  4.12  3.55  25.50  38.53 1.00     3823     3138
##    y0[13]  30.38  30.47  4.01  3.60  23.49  36.64 1.00     4339     3596
##    y0[14]  40.65  41.06  4.19  3.89  33.37  46.98 1.00     2191     3456
##    y0[15]  28.60  28.10  4.52  4.47  22.19  36.61 1.01      664     2621
##    y0[16]  42.88  42.35  4.53  4.39  36.17  50.94 1.01      886     3340
##    y0[17]  15.05  14.72  4.25  3.93   8.64  22.68 1.00     1252     2241
##    y0[18]   8.06   7.66  4.26  3.94   1.48  15.46 1.00     1402     3120
##    y0[19]   5.28   5.49  4.18  3.70  -1.74  11.86 1.00     2682     2887
##    y0[20]  27.62  28.11  4.25  4.09  20.13  33.75 1.01      884     3314
##    m1[1]    4.47   4.42  1.95  1.84   1.32   7.70 1.01     1037     1453
##    m1[2]    8.96   8.96  1.65  1.55   6.29  11.72 1.01     1108     1499
##    m1[3]   13.45  13.45  1.38  1.30  11.24  15.76 1.01     1228     1564
##    m1[4]   17.94  17.96  1.17  1.10  16.07  19.86 1.00     1437     1833
##    m1[5]   22.43  22.46  1.05  1.00  20.75  24.13 1.00     1761     1962
##    m1[6]   26.92  26.95  1.05  0.96  25.23  28.64 1.00     2009     2082
##    m1[7]   31.42  31.44  1.17  1.08  29.54  33.23 1.00     1821     2022
##    m1[8]   35.91  35.96  1.38  1.26  33.72  38.01 1.00     1464     1929
##    m1[9]   40.40  40.46  1.64  1.51  37.74  42.92 1.00     1211     1494
##    m1[10]  44.89  44.96  1.94  1.78  41.69  47.93 1.00     1096     1518
##    x1[1]   -1.18  -1.18  1.64  0.98  -3.91   1.55 1.04     4015     1415
##    x1[2]    1.38   1.41  1.63  1.01  -1.43   4.15 1.04     4114     1970
##    x1[3]    4.05   4.01  1.69  0.99   1.35   6.95 1.03     3549     1609
##    x1[4]    6.61   6.59  1.63  1.00   3.93   9.36 1.03     3677     2227
##    x1[5]    9.15   9.15  1.69  1.05   6.34  11.83 1.04     4111     2304
##    x1[6]   11.70  11.70  1.64  0.96   8.99  14.38 1.03     3896     1756
##    x1[7]   14.32  14.29  1.68  1.05  11.58  17.14 1.04     3890     2313
##    x1[8]   16.92  16.91  1.63  1.00  14.20  19.61 1.04     4201     2090
##    x1[9]   19.49  19.47  1.69  0.99  16.79  22.30 1.04     3808     2045
##    x1[10]  22.02  22.05  1.67  0.99  19.24  24.76 1.04     3835     2267
##    y1[1]    4.48   4.43  4.15  3.74  -2.13  11.30 1.00     2580     2745
##    y1[2]    8.96   8.86  4.02  3.59   2.65  15.46 1.00     2525     3112
##    y1[3]   13.49  13.47  3.84  3.34   7.36  19.88 1.00     3089     3219
##    y1[4]   18.01  18.00  3.82  3.40  11.84  24.38 1.00     3381     3288
##    y1[5]   22.49  22.42  3.76  3.13  16.32  28.69 1.01     3807     3476
##    y1[6]   26.87  26.87  3.72  3.30  20.60  32.98 1.00     3802     2987
##    y1[7]   31.39  31.40  3.79  3.27  25.27  37.70 1.00     3324     3011
##    y1[8]   35.86  35.86  3.87  3.35  29.54  42.21 1.01     3019     2990
##    y1[9]   40.39  40.42  3.97  3.48  33.89  46.91 1.01     3022     2701
##    y1[10]  44.90  44.92  4.11  3.63  38.26  51.76 1.00     2757     2626

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -34076.2 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24      -32.6737   0.000167123   0.000609913      0.8415      0.8415       58    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -32.67
##    b0         6.72
##    b1         1.82
##    s          3.11
##    m0[1]     12.17
##    m0[2]     22.76
##    m0[3]     10.33
##    m0[4]     18.05
##    m0[5]     16.65
##    m0[6]      9.37
##    m0[7]     19.72
##    m0[8]     12.30
##    m0[9]     11.54
##    m0[10]    20.39
##    m0[11]     8.59
##    m0[12]    22.16
##    m0[13]    17.16
##    m0[14]    19.07
##    m0[15]     8.69
##    m0[16]    20.91
##    m0[17]    14.00
##    m0[18]    12.49
##    m0[19]    15.69
##    m0[20]     9.98
##    y0[1]     15.32
##    y0[2]     25.10
##    y0[3]      9.12
##    y0[4]     19.79
##    y0[5]     18.27
##    y0[6]      8.03
##    y0[7]     16.79
##    y0[8]      6.34
##    y0[9]     13.10
##    y0[10]    22.78
##    y0[11]     6.79
##    y0[12]    23.02
##    y0[13]    13.50
##    y0[14]    15.23
##    y0[15]    12.08
##    y0[16]    21.69
##    y0[17]    17.13
##    y0[18]    13.76
##    y0[19]    12.63
##    y0[20]    12.74
##    m1[1]      8.59
##    m1[2]     10.16
##    m1[3]     11.74
##    m1[4]     13.31
##    m1[5]     14.89
##    m1[6]     16.46
##    m1[7]     18.04
##    m1[8]     19.61
##    m1[9]     21.19
##    m1[10]    22.76
##    y1[1]     10.11
##    y1[2]      7.40
##    y1[3]     11.55
##    y1[4]      9.03
##    y1[5]     16.79
##    y1[6]     16.39
##    y1[7]     15.51
##    y1[8]     15.82
##    y1[9]     24.80
##    y1[10]    23.42
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.14 -32.79 1.31 1.07 -35.64 -31.74 1.01      493      637
##    b0       6.66   6.74 1.60 1.57   3.94   9.11 1.01      311      301
##    b1       1.83   1.82 0.30 0.29   1.36   2.32 1.01      406      472
##    s        3.52   3.42 0.66 0.63   2.65   4.75 1.00     1236     1347
##    m0[1]   12.15  12.18 0.95 0.92  10.60  13.64 1.01      389      423
##    m0[2]   22.80  22.81 1.49 1.44  20.34  25.28 1.00     1075     1178
##    m0[3]   10.30  10.35 1.13 1.11   8.38  12.08 1.01      331      346
##    m0[4]   18.07  18.07 0.95 0.89  16.44  19.60 1.00     2043     1281
##    m0[5]   16.65  16.65 0.85 0.80  15.22  18.01 1.00     1644     1248
##    m0[6]    9.33   9.38 1.25 1.21   7.23  11.27 1.01      319      322
##    m0[7]   19.74  19.74 1.11 1.04  17.88  21.56 1.00     1655     1353
##    m0[8]   12.27  12.29 0.94 0.91  10.75  13.78 1.01      395      454
##    m0[9]   11.51  11.55 1.01 0.98   9.85  13.11 1.01      360      381
##    m0[10]  20.42  20.42 1.19 1.14  18.45  22.35 1.00     1462     1247
##    m0[11]   8.54   8.61 1.35 1.32   6.25  10.60 1.01      314      305
##    m0[12]  22.20  22.20 1.41 1.37  19.89  24.52 1.00     1148     1121
##    m0[13]  17.17  17.17 0.88 0.82  15.67  18.57 1.00     1884     1298
##    m0[14]  19.09  19.08 1.04 0.98  17.34  20.79 1.00     1850     1281
##    m0[15]   8.64   8.71 1.33 1.30   6.37  10.68 1.01      314      310
##    m0[16]  20.94  20.94 1.25 1.20  18.87  22.97 1.00     1351     1240
##    m0[17]  13.99  13.99 0.84 0.80  12.63  15.33 1.00      574      886
##    m0[18]  12.47  12.48 0.92 0.91  10.96  13.95 1.01      407      449
##    m0[19]  15.69  15.69 0.82 0.77  14.30  17.00 1.00     1152     1186
##    m0[20]   9.94  10.00 1.17 1.15   7.95  11.77 1.01      326      324
##    y0[1]   12.14  12.18 3.77 3.51   5.96  18.32 1.00     1650     1893
##    y0[2]   22.72  22.82 3.83 3.62  16.35  29.05 1.00     1907     1881
##    y0[3]   10.23  10.21 3.76 3.47   4.11  16.37 1.00     1314     1609
##    y0[4]   18.21  18.20 3.67 3.57  12.41  24.24 1.00     1961     1703
##    y0[5]   16.64  16.61 3.76 3.77  10.73  22.71 1.00     2086     1822
##    y0[6]    9.49   9.51 3.79 3.67   3.22  15.53 1.00     1450     1366
##    y0[7]   19.64  19.63 3.78 3.50  13.51  25.80 1.00     2064     1966
##    y0[8]   12.12  12.10 3.78 3.54   6.11  18.39 1.00     1785     1905
##    y0[9]   11.55  11.47 3.74 3.56   5.66  17.51 1.00     1579     1643
##    y0[10]  20.40  20.38 3.86 3.64  14.17  26.61 1.00     1740     1849
##    y0[11]   8.55   8.54 3.85 3.71   2.14  14.77 1.00     1139     1591
##    y0[12]  22.12  22.08 3.77 3.61  15.94  28.49 1.00     1727     1961
##    y0[13]  17.17  17.20 3.81 3.60  10.52  23.18 1.00     1988     1971
##    y0[14]  19.14  19.13 3.77 3.69  13.03  25.42 1.00     2137     1822
##    y0[15]   8.68   8.73 3.97 3.94   2.02  15.29 1.00     1382     1380
##    y0[16]  20.95  20.90 3.90 3.79  14.60  27.54 1.00     2005     1633
##    y0[17]  14.00  13.85 3.73 3.38   8.00  20.19 1.00     1844     1807
##    y0[18]  12.54  12.44 3.83 3.59   6.35  18.76 1.00     1689     1966
##    y0[19]  15.60  15.66 3.62 3.39   9.70  21.49 1.00     1991     1658
##    y0[20]   9.93  10.02 3.71 3.45   3.69  15.84 1.00     1340     1697
##    m1[1]    8.54   8.61 1.35 1.32   6.25  10.60 1.01      314      305
##    m1[2]   10.13  10.18 1.15 1.13   8.18  11.93 1.01      329      336
##    m1[3]   11.71  11.74 0.99 0.97  10.09  13.28 1.01      368      390
##    m1[4]   13.30  13.30 0.87 0.85  11.87  14.68 1.01      478      613
##    m1[5]   14.88  14.88 0.82 0.78  13.53  16.20 1.00      789     1107
##    m1[6]   16.46  16.46 0.85 0.79  15.04  17.80 1.00     1549     1298
##    m1[7]   18.05  18.05 0.95 0.89  16.43  19.58 1.00     2044     1281
##    m1[8]   19.63  19.63 1.10 1.03  17.79  21.44 1.00     1687     1353
##    m1[9]   21.22  21.21 1.29 1.24  19.11  23.30 1.00     1298     1239
##    m1[10]  22.80  22.81 1.49 1.44  20.34  25.28 1.00     1075     1178
##    y1[1]    8.73   8.82 3.78 3.58   2.54  14.87 1.00     1159     1515
##    y1[2]   10.15  10.17 3.76 3.57   3.95  16.15 1.00     1435     1580
##    y1[3]   11.64  11.59 3.75 3.65   5.63  17.58 1.00     1866     1842
##    y1[4]   13.24  13.19 3.64 3.43   7.33  19.35 1.00     1619     1582
##    y1[5]   14.94  14.93 3.69 3.63   8.94  20.96 1.00     1832     1763
##    y1[6]   16.46  16.35 3.68 3.51  10.59  22.41 1.00     1912     1854
##    y1[7]   17.97  17.95 3.61 3.37  11.94  23.97 1.00     1896     1806
##    y1[8]   19.66  19.67 3.73 3.58  13.51  25.84 1.00     1849     1861
##    y1[9]   21.12  21.06 3.91 3.56  14.67  27.45 1.00     1744     1675
##    y1[10]  22.81  22.82 3.88 3.58  16.53  29.10 1.00     1383     1582

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -100.194 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -13.2521    1.9685e-05    0.00237116      0.3438      0.3438       28    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -13.25
##    b0         5.49
##    b1         1.95
##    s          0.58
##    m0[1]     11.33
##    m0[2]     22.65
##    m0[3]      9.36
##    m0[4]     17.62
##    m0[5]     16.11
##    m0[6]      8.32
##    m0[7]     19.40
##    m0[8]     11.46
##    m0[9]     10.65
##    m0[10]    20.12
##    m0[11]     7.49
##    m0[12]    22.01
##    m0[13]    16.66
##    m0[14]    18.70
##    m0[15]     7.60
##    m0[16]    20.67
##    m0[17]    13.28
##    m0[18]    11.67
##    m0[19]    15.09
##    m0[20]     8.98
##    y0[1]     10.98
##    y0[2]     23.16
##    y0[3]      8.17
##    y0[4]     10.04
##    y0[5]     17.18
##    y0[6]      7.49
##    y0[7]     20.39
##    y0[8]     11.94
##    y0[9]     11.27
##    y0[10]    20.19
##    y0[11]     0.69
##    y0[12]    22.09
##    y0[13]    16.35
##    y0[14]    25.24
##    y0[15]     8.13
##    y0[16]    21.04
##    y0[17]    13.56
##    y0[18]    11.65
##    y0[19]    13.20
##    y0[20]     9.77
##    m1[1]      7.49
##    m1[2]      9.17
##    m1[3]     10.86
##    m1[4]     12.54
##    m1[5]     14.23
##    m1[6]     15.91
##    m1[7]     17.60
##    m1[8]     19.28
##    m1[9]     20.97
##    m1[10]    22.65
##    y1[1]      7.18
##    y1[2]      9.58
##    y1[3]      3.83
##    y1[4]     12.06
##    y1[5]     13.83
##    y1[6]     16.78
##    y1[7]     17.07
##    y1[8]     19.85
##    y1[9]     20.83
##    y1[10]    22.43
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.55 -15.21   1.39 1.16 -18.21 -13.99 1.00      625      862
##    b0       5.46   5.46   0.51 0.49   4.63   6.29 1.01      494      555
##    b1       1.95   1.95   0.09 0.09   1.80   2.09 1.00      592      884
##    s        0.76   0.73   0.25 0.24   0.40   1.21 1.00      598      386
##    m0[1]   11.30  11.31   0.30 0.29  10.79  11.79 1.01      586      675
##    m0[2]   22.63  22.63   0.43 0.37  21.91  23.32 1.00     1168     1134
##    m0[3]    9.33   9.34   0.36 0.33   8.74   9.92 1.01      522      622
##    m0[4]   17.59  17.60   0.28 0.25  17.13  18.03 1.00     1834     1242
##    m0[5]   16.09  16.09   0.25 0.23  15.66  16.48 1.00     1483     1316
##    m0[6]    8.30   8.30   0.40 0.37   7.64   8.93 1.01      506      574
##    m0[7]   19.38  19.38   0.32 0.28  18.84  19.87 1.00     1630     1346
##    m0[8]   11.43  11.44   0.30 0.28  10.93  11.91 1.01      594      675
##    m0[9]   10.62  10.63   0.32 0.30  10.09  11.14 1.01      556      634
##    m0[10]  20.10  20.10   0.34 0.30  19.53  20.62 1.00     1512     1265
##    m0[11]   7.46   7.47   0.43 0.40   6.75   8.16 1.01      499      565
##    m0[12]  21.99  21.99   0.41 0.35  21.31  22.63 1.00     1238     1107
##    m0[13]  16.64  16.64   0.26 0.23  16.19  17.05 1.00     1647     1299
##    m0[14]  18.68  18.68   0.30 0.26  18.18  19.14 1.00     1768     1333
##    m0[15]   7.57   7.57   0.42 0.39   6.86   8.25 1.01      499      561
##    m0[16]  20.65  20.65   0.36 0.31  20.05  21.21 1.00     1432     1195
##    m0[17]  13.25  13.26   0.26 0.24  12.81  13.68 1.01      767      811
##    m0[18]  11.64  11.65   0.29 0.28  11.15  12.12 1.01      608      675
##    m0[19]  15.06  15.06   0.25 0.23  14.63  15.46 1.00     1178     1217
##    m0[20]   8.95   8.96   0.37 0.35   8.33   9.56 1.01      516      635
##    y0[1]   11.29  11.29  15.01 1.20   5.85  15.68 1.00     2057     1942
##    y0[2]   22.45  22.62  21.64 1.14  17.93  27.19 1.00     1936     1883
##    y0[3]   10.11   9.38  75.82 1.19   4.71  14.77 1.00     1873     1965
##    y0[4]   17.72  17.55  17.06 1.21  11.91  22.59 1.00     1878     1802
##    y0[5]   17.10  16.11  25.82 1.13  11.40  21.05 1.00     1809     1753
##    y0[6]    8.57   8.29  16.23 1.16   3.45  12.97 1.00     1933     1655
##    y0[7]   20.95  19.35  47.55 1.15  14.78  24.41 1.00     1748     1545
##    y0[8]    7.32  11.45 157.92 1.19   6.45  16.60 1.00     1915     1894
##    y0[9]    7.97  10.61 303.88 1.22   5.28  15.64 1.00     1890     1799
##    y0[10]  19.99  20.03  23.97 1.17  14.84  24.56 1.00     1715     1774
##    y0[11]   7.41   7.51  12.27 1.22   3.15  12.09 1.00     1769     2097
##    y0[12]  22.19  21.96  20.59 1.15  17.69  26.58 1.00     1991     1929
##    y0[13]  24.35  16.64 326.61 1.13  12.65  21.32 1.00     2140     1933
##    y0[14]  18.95  18.68  27.92 1.14  13.77  23.91 1.00     2059     1894
##    y0[15]   7.90   7.54  79.37 1.24   1.85  12.27 1.00     1636     1905
##    y0[16]  21.45  20.66  17.26 1.13  16.37  24.86 1.00     1874     1823
##    y0[17]  13.66  13.22  35.95 1.19   7.68  17.98 1.00     1807     1643
##    y0[18]  11.21  11.59  20.90 1.11   6.38  15.87 1.00     2180     2037
##    y0[19]  15.40  15.06  22.02 1.14  10.38  20.18 1.00     1978     1865
##    y0[20]   7.02   8.95  78.36 1.18   3.45  13.90 1.00     1907     1894
##    m1[1]    7.46   7.47   0.43 0.40   6.75   8.16 1.01      499      565
##    m1[2]    9.15   9.15   0.37 0.34   8.55   9.74 1.01      519      635
##    m1[3]   10.83  10.84   0.31 0.30  10.31  11.34 1.01      564      647
##    m1[4]   12.52  12.52   0.27 0.26  12.06  12.96 1.01      681      745
##    m1[5]   14.20  14.21   0.25 0.24  13.78  14.62 1.00      942      954
##    m1[6]   15.89  15.89   0.25 0.23  15.47  16.28 1.00     1433     1256
##    m1[7]   17.57  17.58   0.27 0.25  17.11  18.01 1.00     1833     1242
##    m1[8]   19.26  19.26   0.32 0.27  18.73  19.75 1.00     1653     1346
##    m1[9]   20.94  20.95   0.37 0.31  20.33  21.53 1.00     1388     1135
##    m1[10]  22.63  22.63   0.43 0.37  21.91  23.32 1.00     1168     1134
##    y1[1]    5.10   7.47 163.26 1.24   3.06  12.61 1.01     1460     1997
##    y1[2]    7.92   9.13  51.20 1.14   5.34  13.25 1.00     1651     1835
##    y1[3]   10.32  10.84  12.90 1.09   6.68  15.52 1.00     1999     1908
##    y1[4]   12.98  12.53  25.21 1.13   7.52  17.34 1.00     1974     1846
##    y1[5]   13.85  14.20  10.99 1.11   9.00  18.71 1.00     2036     1959
##    y1[6]   16.96  15.89  42.38 1.16  12.01  20.33 1.00     1904     1974
##    y1[7]   17.51  17.58  83.77 1.17  12.98  21.73 1.00     1925     1792
##    y1[8]   19.24  19.27  16.73 1.08  14.91  23.52 1.00     1922     1671
##    y1[9]   14.26  20.95 328.50 1.20  16.27  25.62 1.00     1944     1845
##    y1[10]  23.20  22.64  24.25 1.25  17.89  28.18 1.00     2077     1872

censored data

objective variable has NA

(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)

ex11-0.stan

data{
  int N0;
  array[N0] vector[2] xy;
  int N1;
  vector[N1] x1;
}
parameters{
  vector[2] m;
  cov_matrix[2] s;
}
model{
  target+=multi_normal_lpdf(xy | m, s);
  x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
  vector[2] xy1;
  xy1=multi_normal_rng(m,s);
  real cr;
  cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}
n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)
## [1] 0.889
qplot(x,y)

L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)
## [1] 0.789
qplot(x0,y0)

mdl=cmdstan_model('./ex11-0.stan') 

data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -1794.89 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       49        -127.4   0.000806172    0.00104209      0.8332      0.8332       70    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    -127.40
##    m[1]       5.73
##    m[2]      27.80
##    s[1,1]     5.00
##    s[2,1]    17.09
##    s[1,2]    17.09
##    s[2,2]    73.08
##    xy1[1]     6.13
##    xy1[2]    27.00
##    cr         0.89
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.5 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -123.59 -123.24  1.75  1.52 -126.99 -121.47 1.01      543      652
##    m[1]      5.74    5.74  0.46  0.47    5.00    6.50 1.01      573      639
##    m[2]     27.84   27.87  1.89  1.93   24.75   30.85 1.01      506      599
##    s[1,1]    6.63    6.24  2.11  1.74    4.06   10.40 1.01      858      840
##    s[2,1]   22.89   21.49  8.70  7.19   12.06   38.83 1.01      526      516
##    s[1,2]   22.89   21.49  8.70  7.19   12.06   38.83 1.01      526      516
##    s[2,2]  101.75   94.20 41.53 33.51   51.12  179.53 1.02      523      523
##    xy1[1]    5.82    5.76  2.51  2.47    1.65    9.96 1.00     1864     1745
##    xy1[2]   28.02   27.99  9.82  9.31   11.46   44.56 1.00     1858     1758
##    cr        0.88    0.89  0.07  0.05    0.76    0.95 1.00      617      866

xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])
## [1] 0.871
qplot(xy[,,1],xy[,,2])

objective variable is censored

y i=1-N, y~N(m,s)  
  actual          ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
mdl=cmdstan_model('./ex8-0.stan')

data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -566.386 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       29      -10.3006   0.000345667   0.000330398      0.9893      0.9893       41    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -10.30
##    b0        13.44
##    b1         2.02
##    s          1.55
##    m0[1]     25.98
##    m0[2]     23.61
##    m0[3]     21.55
##    m0[4]     19.65
##    m0[5]     16.19
##    m0[6]     21.79
##    m0[7]     23.68
##    m0[8]     16.39
##    m0[9]     28.83
##    m0[10]    26.16
##    m0[11]    26.50
##    y0[1]     23.60
##    y0[2]     21.30
##    y0[3]     22.02
##    y0[4]     22.74
##    y0[5]     13.60
##    y0[6]     21.89
##    y0[7]     23.99
##    y0[8]     17.77
##    y0[9]     29.19
##    y0[10]    24.87
##    y0[11]    25.03
##    m1[1]     13.55
##    m1[2]     15.47
##    m1[3]     17.40
##    m1[4]     19.33
##    m1[5]     21.26
##    m1[6]     23.18
##    m1[7]     25.11
##    m1[8]     27.04
##    m1[9]     28.97
##    m1[10]    30.89
##    y1[1]     12.20
##    y1[2]     16.28
##    y1[3]     18.75
##    y1[4]     18.52
##    y1[5]     19.20
##    y1[6]     21.97
##    y1[7]     26.35
##    y1[8]     27.96
##    y1[9]     29.76
##    y1[10]    33.49
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -11.70 -11.29 1.56 1.23 -14.76 -10.05 1.01      529      632
##    b0      13.27  13.30 1.70 1.57  10.45  15.86 1.02      319      318
##    b1       2.06   2.05 0.34 0.32   1.53   2.62 1.01      361      388
##    s        2.03   1.90 0.64 0.48   1.30   3.29 1.00      726      905
##    m0[1]   26.03  26.00 0.86 0.72  24.68  27.37 1.00     1448     1373
##    m0[2]   23.62  23.59 0.68 0.58  22.58  24.70 1.00     1683     1288
##    m0[3]   21.52  21.52 0.70 0.64  20.44  22.61 1.00      663      846
##    m0[4]   19.59  19.60 0.84 0.76  18.22  20.92 1.01      405      491
##    m0[5]   16.06  16.11 1.29 1.18  13.97  18.06 1.02      326      328
##    m0[6]   21.76  21.76 0.69 0.62  20.69  22.84 1.00      738      872
##    m0[7]   23.68  23.65 0.68 0.58  22.64  24.77 1.00     1712     1272
##    m0[8]   16.27  16.31 1.26 1.16  14.21  18.22 1.02      327      328
##    m0[9]   28.92  28.90 1.22 1.02  26.99  30.95 1.00      733      894
##    m0[10]  26.21  26.18 0.88 0.74  24.85  27.59 1.00     1369     1373
##    m0[11]  26.55  26.52 0.92 0.78  25.14  27.99 1.00     1258     1246
##    y0[1]   26.01  26.02 2.33 2.11  22.14  29.66 1.00     1879     1933
##    y0[2]   23.62  23.66 2.28 2.03  19.96  27.17 1.00     1699     1843
##    y0[3]   21.54  21.55 2.25 2.07  17.90  25.16 1.00     1402     1641
##    y0[4]   19.45  19.52 2.34 2.11  15.78  23.05 1.00     1007     1165
##    y0[5]   16.04  16.13 2.49 2.30  11.95  19.97 1.00      765      925
##    y0[6]   21.66  21.69 2.28 1.97  17.95  25.15 1.00     1560     1595
##    y0[7]   23.66  23.67 2.25 1.98  20.12  27.27 1.00     1847     1916
##    y0[8]   16.31  16.38 2.46 2.22  12.23  20.31 1.00      769     1233
##    y0[9]   28.88  28.83 2.40 2.16  25.14  32.80 1.00     1453     1571
##    y0[10]  26.20  26.18 2.28 2.01  22.49  29.89 1.00     1745     1767
##    y0[11]  26.58  26.47 2.34 2.10  22.95  30.40 1.00     1695     1778
##    m1[1]   13.38  13.41 1.68 1.55  10.59  15.95 1.02      319      318
##    m1[2]   15.34  15.39 1.39 1.29  13.06  17.48 1.02      323      322
##    m1[3]   17.30  17.32 1.12 1.02  15.46  19.04 1.02      336      338
##    m1[4]   19.26  19.28 0.88 0.80  17.82  20.65 1.01      387      468
##    m1[5]   21.22  21.23 0.71 0.64  20.11  22.33 1.00      592      702
##    m1[6]   23.18  23.15 0.67 0.57  22.16  24.24 1.00     1455     1271
##    m1[7]   25.14  25.11 0.78 0.66  23.93  26.35 1.00     1755     1379
##    m1[8]   27.10  27.07 0.98 0.83  25.58  28.64 1.00     1060     1143
##    m1[9]   29.06  29.04 1.24 1.04  27.09  31.14 1.00      719      885
##    m1[10]  31.02  31.01 1.52 1.33  28.60  33.51 1.00      591      740
##    y1[1]   13.42  13.44 2.74 2.44   8.82  17.67 1.01      562      766
##    y1[2]   15.31  15.37 2.61 2.37  11.08  19.36 1.00      862     1079
##    y1[3]   17.29  17.30 2.42 2.14  13.43  21.06 1.00      992     1316
##    y1[4]   19.26  19.33 2.31 2.09  15.50  22.88 1.00     1349     1469
##    y1[5]   21.21  21.21 2.29 2.04  17.58  25.02 1.00     1822     1652
##    y1[6]   23.17  23.11 2.24 2.02  19.58  26.73 1.00     1748     1758
##    y1[7]   25.11  25.00 2.23 2.02  21.61  28.61 1.00     1903     1688
##    y1[8]   27.09  27.07 2.30 2.08  23.44  30.77 1.00     1687     1605
##    y1[9]   29.07  28.94 2.48 2.28  25.17  33.19 1.00     1304     1257
##    y1[10]  31.01  30.95 2.62 2.44  26.94  35.23 1.00      997     1316

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -52.8227 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23      -23.4739    0.00204867   0.000134233           1           1       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -23.47
##    b0         8.42
##    b1         2.97
##    s          3.14
##    m0[1]     26.84
##    m0[2]     23.36
##    m0[3]     20.33
##    m0[4]     17.54
##    m0[5]     12.45
##    m0[6]     20.69
##    m0[7]     23.46
##    m0[8]     12.75
##    m0[9]     31.02
##    m0[10]    27.10
##    m0[11]    27.60
##    y0[1]     27.44
##    y0[2]     20.73
##    y0[3]     18.51
##    y0[4]     12.20
##    y0[5]     16.59
##    y0[6]     24.23
##    y0[7]     24.74
##    y0[8]     17.58
##    y0[9]     31.11
##    y0[10]    29.86
##    y0[11]    25.67
##    m1[1]      8.57
##    m1[2]     11.40
##    m1[3]     14.23
##    m1[4]     17.07
##    m1[5]     19.90
##    m1[6]     22.73
##    m1[7]     25.56
##    m1[8]     28.39
##    m1[9]     31.23
##    m1[10]    34.06
##    y1[1]      9.63
##    y1[2]     15.14
##    y1[3]     16.92
##    y1[4]     15.03
##    y1[5]     22.33
##    y1[6]     27.32
##    y1[7]     24.73
##    y1[8]     30.02
##    y1[9]     33.52
##    y1[10]    32.69
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -24.01 -23.62 1.39 1.07 -26.58 -22.52 1.01      581      675
##    b0       7.14   7.37 2.89 2.45   2.05  11.53 1.02      234      366
##    b1       3.23   3.19 0.58 0.51   2.36   4.24 1.01      288      357
##    s        4.16   3.93 1.27 0.98   2.66   6.56 1.01      718      751
##    m0[1]   27.13  27.04 1.44 1.35  25.02  29.53 1.00     1047      804
##    m0[2]   23.35  23.33 1.13 1.06  21.64  25.21 1.00     1387     1059
##    m0[3]   20.07  20.14 1.15 1.04  18.15  21.83 1.01      534      724
##    m0[4]   17.04  17.17 1.41 1.22  14.63  19.18 1.02      280      474
##    m0[5]   11.51  11.71 2.18 1.84   7.77  14.86 1.02      222      360
##    m0[6]   20.45  20.50 1.13 1.02  18.57  22.20 1.01      592      730
##    m0[7]   23.46  23.43 1.13 1.07  21.75  25.33 1.00     1409     1096
##    m0[8]   11.84  12.03 2.13 1.80   8.20  15.12 1.02      221      360
##    m0[9]   31.67  31.52 2.06 1.88  28.65  35.11 1.00      617      618
##    m0[10]  27.42  27.31 1.47 1.39  25.26  29.87 1.00      999      800
##    m0[11]  27.96  27.85 1.54 1.43  25.71  30.54 1.00      917      800
##    y0[1]   27.25  27.15 4.53 4.14  20.03  34.58 1.00     1953     1723
##    y0[2]   23.15  23.22 4.47 3.97  15.96  30.09 1.00     2060     1726
##    y0[3]   20.23  20.18 4.56 3.94  13.01  27.65 1.00     1891     1605
##    y0[4]   17.06  17.26 4.58 3.93   9.15  24.12 1.00     1465     1417
##    y0[5]   11.44  11.68 4.82 4.49   3.53  18.73 1.01      910     1239
##    y0[6]   20.49  20.62 4.52 3.96  12.96  27.51 1.00     2020     1775
##    y0[7]   23.47  23.54 4.49 4.08  16.03  30.55 1.00     1931     1457
##    y0[8]   11.93  12.15 4.86 4.26   3.50  19.56 1.00     1053     1224
##    y0[9]   31.80  31.56 4.87 4.40  24.22  39.83 1.00     1417     1480
##    y0[10]  27.30  27.15 4.68 4.25  19.92  35.02 1.00     1884     1950
##    y0[11]  27.88  27.83 4.77 4.33  20.48  35.38 1.00     1939     1635
##    m1[1]    7.30   7.54 2.86 2.43   2.26  11.66 1.02      233      365
##    m1[2]   10.38  10.60 2.36 1.99   6.26  14.00 1.02      224      354
##    m1[3]   13.45  13.61 1.89 1.61  10.17  16.32 1.02      222      389
##    m1[4]   16.52  16.65 1.47 1.28  13.99  18.74 1.02      263      449
##    m1[5]   19.60  19.68 1.18 1.07  17.63  21.40 1.01      474      662
##    m1[6]   22.67  22.66 1.11 1.01  20.93  24.50 1.00     1200     1002
##    m1[7]   25.74  25.68 1.29 1.24  23.77  27.90 1.00     1290     1070
##    m1[8]   28.82  28.70 1.65 1.53  26.42  31.59 1.00      811      736
##    m1[9]   31.89  31.75 2.10 1.89  28.86  35.39 1.00      606      624
##    m1[10]  34.96  34.75 2.59 2.29  31.20  39.27 1.01      504      502
##    y1[1]    7.27   7.51 5.21 4.73  -1.41  14.91 1.00      792      994
##    y1[2]   10.27  10.45 5.02 4.41   2.07  17.79 1.01      945     1168
##    y1[3]   13.32  13.53 4.71 4.29   5.51  20.61 1.00     1032     1147
##    y1[4]   16.47  16.62 4.65 4.11   8.66  23.56 1.00     1217     1492
##    y1[5]   19.63  19.58 4.43 4.06  12.46  26.72 1.00     1611     1795
##    y1[6]   22.68  22.62 4.60 4.14  15.49  30.21 1.00     1924     1933
##    y1[7]   25.71  25.72 4.53 4.09  18.42  33.08 1.00     1792     1752
##    y1[8]   28.79  28.64 4.67 4.43  21.54  36.39 1.00     1486     1585
##    y1[9]   31.97  31.62 4.82 4.37  24.78  40.29 1.00     1767     1903
##    y1[10]  34.96  34.76 5.12 4.51  27.01  43.68 1.00     1065     1042

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
data=list(N=n,
          x=sample(0:1,n,replace=T),
          y=sample(0:1,n,replace=T))

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -12.8506 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -11.2128   0.000179292    5.9677e-06           1           1        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -11.21
##      p        0.45
##      se       0.78
##      sp       0.73
##      ppv      0.70
##      npv      0.80
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -17.54 -17.20 1.33 1.15 -19.99 -16.05 1.00      767      828
##      p      0.51   0.51 0.30 0.38   0.05   0.96 1.00      460      495
##      se     0.73   0.75 0.13 0.12   0.49   0.92 1.00     1918     1304
##      sp     0.69   0.71 0.12 0.13   0.47   0.87 1.00     3189     1493
##      ppv    0.65   0.73 0.28 0.30   0.10   0.98 1.00      484      570
##      npv    0.65   0.73 0.29 0.31   0.10   0.99 1.00      478      617

bimodal distribution, mixed normal distribution

n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()

EM algorithm

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -254 100  6 -535 -538
## 
## Clustering table:
##  1  2  3 
## 32 38 30
rst$parameters
## $pro
## [1] 0.320 0.384 0.295
## 
## $mean
##       1       2       3 
## -5.1792  0.0431  4.8765 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 1.14
## 
## 
## $Vinv
## NULL
plot(rst)

ex16-1.stan

data {
  int N;
  vector[N] y;;
}
parameters {
  simplex[3] h; //ratio of mix
  real m1;
  real m2;
  real m3;
  real<lower=0> s1;
  real<lower=0> s2;
  real<lower=0> s3;
}
model {
  s1~cauchy(0,5);
  s2~cauchy(0,5);
  s3~cauchy(0,5);
  for (i in 1:N) {
    vector[3] p;
    p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
    p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
    p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
    target+=log_sum_exp(p);
  }
}
mdl=cmdstan_model('./ex16-1.stan')

data=list(N=n,y=y)
mle=mdl$optimize(data=data)
## Initial log joint probability = -635.698 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -265.411   0.000104163     0.0321676      0.3444     0.03444      130    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      100      -265.411   0.000257792     0.0115348           1           1      131    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -265.41
##      h[1]     0.70
##      h[2]     0.25
##      h[3]     0.05
##      m1      -2.18
##      m2       5.05
##      m3       0.71
##      s1       3.10
##      s2       0.88
##      s3       0.04
mcmc=goMCMC(mdl,data,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 1 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 2 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 3 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 1 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 2 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 1 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 1 finished in 1.3 seconds.
## Chain 2 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 3 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 2 finished in 1.3 seconds.
## Chain 3 finished in 1.3 seconds.
## Chain 4 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 4 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 4 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 4 finished in 2.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.5 seconds.
## Total execution time: 2.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -261.34 -260.89 2.58 2.08 -265.73 -258.30 1.00     2237     1424
##      h[1]    0.34    0.34 0.07 0.07    0.24    0.46 1.25       11       94
##      h[2]    0.32    0.32 0.07 0.06    0.23    0.44 1.24       11       40
##      h[3]    0.33    0.33 0.07 0.06    0.24    0.44 1.17       15       72
##      m1     -0.06    0.08 3.71 6.19   -5.34    5.09 2.33        4       26
##      m2      1.14    2.49 4.28 3.71   -5.35    5.18 2.10        5       29
##      m3     -1.34   -2.31 4.16 4.21   -5.47    5.07 2.10        5       26
##      s1      1.17    1.12 0.66 0.22    0.82    1.62 1.05       55     3059
##      s2      1.15    1.10 0.40 0.22    0.81    1.61 1.06       42     1540
##      s3      1.19    1.15 0.29 0.22    0.84    1.67 1.06       39      195